$title  CF_CO2_NH.GMS

* ###################################

* Per capita income, consumption patterns, and CO2 emissions

* Journal of the Association of Environmental and Resource Economists

* Justin Caron and Thibault Fally

* April 2021

*###################################

* SIMULATION FILE with CRIE preferences

*###################################


$set SLASH \
$if %system.filesys% == UNIX $set SLASH /


* ------------------------------------
* 
* code is set up to feed in production/factor/demand shares "fitted" to the NH demand vector

* declare set of demand estimates to use : theta4, notc, tc
$if not set demandest $set demandest tc
* need to choose H (homothetic) or NH here
$if not set spec $set spec NH

* declare supply side specification: noresprod, resprod, resprod_06, resprod_15,  resprod_50, resprod_06_gas15
* affects supply elasticity among other things

* NOTE: RESPROD_06 USES A SUPPLY ELASTICITY OF 0.75
$if not set prodspec $set  prodspec resprod_06
*$if not set prodspec $set  prodspec resprod_15
*$if not set prodspec $set  prodspec resprod_50
*$if not set prodspec $set  prodspec noresprod

$if not set ds $set ds gtap8gas


parameters income_shock, tcost_shock, income_shock_toresourcefactor;
tcost_shock = 1;

* shock income by 1 pct
income_shock = 1.01;


* set this to 1 if productivity of resource factor is to be shocked as well -- 0 if off
income_shock_toresourcefactor = 1;
$if "%prodspec%"=="noresprod"  income_shock_toresourcefactor = 0;


* CHOOSE NB OF ITERATIONS:  (use 20 for final results, though converges pretty well after ~5)
Sets
T iteration  /T1*T20/;

Sets
r country  /
AUS,NZL,CHN,HKG,JPN,KOR,MNG,TWN,KHM,IDN,LAO,MYS,PHL,SGP,THA,VNM,BGD,IND,NPL,PAK,LKA,CAN,
USA,MEX,ARG,BOL,BRA,CHL,COL,ECU,PRY,PER,URY,VEN,CRI,GTM,HND,NIC,PAN,SLV,AUT,BEL,CYP,CZE,
DNK,EST,FIN,FRA,DEU,GRC,HUN,IRL,ITA,LVA,LTU,LUX,MLT,NLD,POL,PRT,SVK,SVN,ESP,SWE,GBR,CHE,
NOR,ALB,BGR,BLR,HRV,ROU,RUS,UKR,KAZ,KGZ,ARM,AZE,GEO,BHR,IRN,KWT,QAT,SAU,TUR,ARE,EGY,MAR,
TUN,CMR,CIV,GHA,NGA,SEN,ETH,KEN,MDG,MWI,MUS,MOZ,TZA,UGA,ZMB,ZWE,BWA,NAM,ZAF,ISR,OMN
/;


alias (r,s);
alias (r,rloop);

Sets
i industry  /
isr,obs,ros,osg,pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,frs,fsh,coa,oil,gas,omn,cmt,omt,vol,mil,pcr,
sgr,ofd,b_t,tex,wap,lea,lum,ppp,p_c,crp,nmm,i_s,nfm,fmp,mvh,otn,ele,ome,omf,ely,wtr,cns,trd,otp,wtp,atp,cmn,ofi
/;

alias (i,j);

Sets
f factor /Land,UnskLab,SkLab,Capital,NatlRes/;


* Loading data:
parameter fittedexp, coeffs, scale, sigma, theta, beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, incelast;


* Load demand estimates from DEMAND.GMS
$gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_rall.gdx

$load  coeffs scale= sigma_scale fittedexp

* And everything else from DATAPREPARATION.GMS

$gdxin estimates%SLASH%co2_data_%demandest%_avg.gdx
$load  beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, incelast


theta(i) = 4;

* -------
* Theta 4:  H, sigma = scale, NH, sigma = sigma * scale
*  here, the average sigma level is recoverd from Mu

* TC: H, back out thetas assuming the average level of theta = 4?

$if "%demandest%"=="tc" theta(i)$coeffs("backed out theta","coeff",i) = coeffs("backed out theta","coeff",i);

* IF USING THETA = 4 IN THE SIMULATIONS
*$if "%demandest%"=="tc" theta(i) = 4;

*$if "%demandest%"=="homotc" theta(i)$coeffs("backed out theta","coeff",i) = coeffs("backed out theta","coeff",i);

*normalize theta's such that the average equals = 4 and min = 1:
theta(i)$(theta(i) = 0)  = 0;
theta(i)$(theta(i) < 0)  = 0;

* before re-norm:
* some theta's are very large, set max =30
theta(i)$(theta(i) > 30)  = 30;


theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));
theta(i)$(theta(i) = 1)  = 1;
theta(i)$(theta(i) < 1)  = 1;
theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));

*  set max =20 for theta 
theta(i)$(theta(i) > 20)  = 20;
* check: does not affect any sectors in bmk spec

sigma(i) =  coeffs("sigma","coeff",i);

parameter  sigma_rescale;
sigma_rescale = (sum(j$sigma(j), 1) + 4 * sum(j, coeffs("mu","coeff",j) )) / sum(j$sigma(j), sigma(j));
display sigma_rescale;

* TC
* rescale sigmas such that average price elasticity (1-sigma) = average mu * theta
$if "%demandest%"=="tc" $if "%spec%"=="NH" sigma(i) = sigma_rescale * coeffs("sigma","coeff",i)

* For H, use the same average sigma (could also recompute using estimated mu's with H prefs..)
$if "%demandest%"=="tc" $if "%spec%"=="H"  sigma(i) = sigma_rescale;

* Theta 4 / 8: 
$if "%demandest%"=="theta4" $if "%spec%"=="NH"  sigma(i) = scale("non-homoth") * coeffs("sigma","coeff",i);
$if "%demandest%"=="theta8" $if "%spec%"=="NH"  sigma(i) = scale("non-homoth") * coeffs("sigma","coeff",i);


$if "%demandest%"=="theta4" $if "%spec%"=="H"  sigma(i) = scale("homoth") ;
$if "%demandest%"=="theta8" $if "%spec%"=="H"  sigma(i) = scale("homoth") ;

$if "%demandest%"=="theta8" theta(i) = 8;

* sectors with no FD which we have declared "intermediate only":
* setting sigma = average sigma.  those are sectors in the not_fd set

set not_fd(i)  sectors with no sigma estimate (where dropped from estimation);
not_fd(i)$(not sigma(i))  = yes;

* if no sigma, assign average sigma
sigma(i)$not_fd(i)  = sum(i.local$sigma(i), sigma(i)) / sum(i.local$sigma(i), 1) ;

* set theta = 4 for sectors for which we have no estimates
theta(i)$(not theta(i)) = 4;

* ----
* CHOICE OF PRODUCTION SHARES 

* compute fitted production shares from fitted demand shares
* keep observed values for all non-demand related values

parameter finaldemand;

finaldemand("true",i,r) = fd(i,r); 
finaldemand("non-homoth",i,r) = fittedexp("non-homoth",i,r); 
finaldemand("homoth",i,r) = fittedexp("homoth",i,r); 
 

** redefining absorption and production for consistency:
parameters production, absorption;

production(i,r) = sum(s, trade("data",i,r,s));
* INSTEAD OF:	 
*production(i,r) = vom(i,r);
		 
absorption(j,r) = fd(j,r) + sum(i,intersh(j,i,r) * production(i,r));
* INSTEAD OF:
*absorption(i,r) = sum(j, vdfm(i,j,r)+ vifm(i,j,r));
* THIS GUARANTEES WELL DEFINED GOSH COEFFICIENTS -- BUT CONSISTENCY ISSUES WITH IMPORT FLOWS.


display absorption, production;
* Production: some zeros
* Absorption: zeros with previous definition


parameter
dem_sh(i,r)  share of sector i in total final demand
trade_sh(i,r,s)  share of country r in country s's absorption (Import shares)
prod_sh(i,r,s)  share of country s in country r's sales
fact_sh(f,i,r)   share of sector i in demand for factor f
endow_sh(f,r) share of factor f in total endowment (income)
gosh_m(j,i,r) gosh matrix coefficient (how much j is used in industry i)
gosh_f(j,r) gosh matrix coefficient (how much j is used by final consumers);


parameter prod_sh_adj;
parameter vom_adj , vfm_adj;


set fit /non-homoth, homoth, true/;
parameter fitted_prod_sh(i,r,s,fit), fitted_production_absorption;
$gdxin estimates%SLASH%fitted_prod_shares_%spec%_%demandest%.gdx
$load  fitted_prod_sh fitted_production_absorption

* note: HOMOTHETIC not done -- would be easy
*  also wand to adjust demand shares, employement shares, and gosh coeffs

** SPECIFICATION, chose here:
* 1 - bilateral production shares
*prod_sh_adj(i,r,s) = fitted_prod_sh(i,r,s,"true");
prod_sh_adj(i,r,s) = fitted_prod_sh(i,r,s,"non-homoth");
* tested that using only fitted production shares works: eq solutions is 1

* 2 - demand shares
dem_sh(i,r) =  finaldemand("non-homoth",i,r)  / sum(i.local, finaldemand("non-homoth",i,r) );


* 3- absorption and production  --> WILL BE USED TO COMPUTE GOSH COEFFICIENTS 
* keeping observed:
*production(i,r) = production(i,r);
*absorption(j,r) = absorption(j,r);

* or fitted 
production(i,r) = fitted_production_absorption("production","non-homoth",i,r);
absorption(j,r) = fitted_production_absorption("absorption","non-homoth",j,r);

* 4- FACTOR SHARES: also require production and absorption to be computed
* keeping observed (use this to minimize distortions)
*vom_adj(i,r) = vom(i,r);
*vfm_adj(f,i,r) = vfm(f,i,r);

* using 'fitted' absorption/production 
vom_adj(i,r) =  fitted_production_absorption("production","non-homoth",i,r);
vfm_adj(f,i,r) = beta(f,i, r) * vom_adj(i,r);

* -------------------------------------------------------------------------


parameters CHECK_prod_sh;
CHECK_prod_sh(i,r,fit) = round(sum(s,fitted_prod_sh(i,r,s,fit))-1,13);

display CHECK_prod_sh;

* j being used in industry i:
gosh_m(j,i,r)$absorption(j,r) = intersh(j,i,r) * production(i,r) / absorption(j,r);

gosh_f(j,r) = 1 - sum(i,gosh_m(j,i,r));
* normalized to unity when absorption is zero

* check that there is no negative value:
display gosh_f;
* Some negative values because of rounding errors, not a problem (?)

* defined as: from r to s:
trade_sh(i,r,s)$sum(r.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(r.local , trade("data",i,r,s));

parameter comp_ab1, comp_ab2;
comp_ab1(i,s)$absorption(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
comp_ab2(i,s)$(round(sum(r,trade_sh(i,r,s)) - 1,14)) = absorption(i,s);
display comp_ab1, comp_ab2;
* dataset not quite consistent but very close

trade_sh(i,s,s)$(sum(r.local,trade("data",i,r,s))=0)  = 1;

* this is to be replaced if using fitted shares
* defined as: from r to s:
prod_sh(i,r,s)$sum(s.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(s.local , trade("data",i,r,s));

* ADJUSTMENT WHEN SUM IS ZERO:
prod_sh(i,r,r)$(sum(s.local , trade("data",i,r,s))=0)  = 1;

fact_sh(f,i,r) =  beta(f,i, r) * vom_adj(i,r) / sum(i.local,  beta(f,i,r) * vom_adj(i,r))  ;

endow_sh(f,r) = sum(i,vfm_adj(f,i,r)) / sum((f.local,i),vfm_adj(f,i,r)) ;


* CHECKS TO ENSURE THAT CHG.L=1 IS A SOLUTION:

parameters CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;
CHECK_budget(r)= round(sum(i,dem_sh(i,r)) - 1,14);
CHECK_absorption(j,r) = round(gosh_f(j,r) + sum(i$gosh_m(j,i,r),gosh_m(j,i,r)) - 1,14);
CHECK_production(i,r) = round(sum(s,prod_sh(i,r,s)) - 1,14);
CHECK_supplier(i,r) = round(prod(f$beta(f,i,r),1**beta(f,i,r))**theta(i) * prod(j,1**intersh(j,i,r))**theta(i) - 1,14);
CHECK_index(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
CHECK_fmc(f,r) = round(sum(i,fact_sh(f,i,r)) - 1,14);
CHECK_income(r) = round(sum(f,endow_sh(f,r)) - 1,14);

display CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;
* All fine 

parameter trade_sh_adj;
trade_sh_adj(i,r,s)$absorption(i,s) = trade("data",i,r,s)/absorption(i,s);

*---------------------------------------------------------------------------------------------------

* ADJUSTMENT TO GET A CES PROD FUNCTION IN NATURAL RESOURCES:

Set	f_noRes(f)	 factor other THAN NatlRes
	res_sect(i)   sector with sector specfic resource 
	pe(i)		 primary energy 
		;

f_noRes(f) = yes;
f_noRes("NatlRes")  = no;
pe("gas") = yes;
pe("oil") = yes;
pe("coa") = yes;

* ADJUSTMENTS 1 -- for sectors to take out the natural resource factor for non-enegy sectors INCLUDES OMN FSH FRS
* transfer that resource to capital
* i.e, only keep NATLRES for the ENE sectors

set nonatlres(i) / OMN, FSH, FRS /;

* change vfm and beta
vfm_adj("Capital",nonatlres,r) = vfm_Adj("Capital",nonatlres,r) + vfm_adj("NatlRes",nonatlres,r) ;
beta("Capital",nonatlres,r) = beta("Capital",nonatlres,r) + beta("NatlRes",nonatlres,r) ;
vfm_adj("NatlRes",nonatlres,r) = 0;
 beta("NatlRes",nonatlres,r) = 0;

* re-compute factor and endowment shares
fact_sh(f,i,r)$sum(i.local,  beta(f,i,r) * vom_adj(i,r))  =  beta(f,i, r) * vom_adj(i,r) / sum(i.local,  beta(f,i,r) * vom_adj(i,r))  ;
endow_sh(f,r) = sum(i,vfm_adj(f,i,r)) / sum((f.local,i),vfm_adj(f,i,r)) ;


* ADJUSTMENTS 2 -- some countries do not produce certain fossil fuels at all.
* in this case, best to fix W_res to 1
set no_ene_prod(i,r);
no_ene_prod(pe,r)$(not production(pe,r)) = yes;

parameter beta_all(i,r), beta_Res(i,r), beta_fnores(f,i,r), endow_sh_res  ;

* note: as of now, the beta are NOT COUNTRY-SPECIFIC

beta_all(i,r) = sum(f, beta(f,i,r)) ;

beta_fnores(f,i,r)$sum(f.local$f_noRes(f), beta(f,i,r)) = beta(f,i,r) / sum(f.local$f_noRes(f), beta(f,i,r));

beta_Res(i,r)$sum(f.local, beta(f,i,r)) = beta("NatlRes",i,r) / sum(f.local, beta(f,i,r));


parameter	supply_elast(i,r)    targeted elasticity of supply
		nu(i,r)  elasticity of subsitution between resource and non-resource factor  -- to be calibrated to reach desired  supply elasticity ;


supply_elast("oil",r) = 0.75 ;

* necessary to keep nu positive - cannot go much below this
supply_elast("gas",r) = 0.75;
supply_elast("coa",r) = 0.8;

* minimum value for coal to have nu>0 is 0.8. set other two to 0.75
$if "%prodspec%"=="resprod_06" supply_elast(pe,r) = 0.75 ; supply_elast("coa",r) = 0.8;
$if "%prodspec%"=="resprod_15" supply_elast(pe,r) = 1.5;
$if "%prodspec%"=="resprod_50" supply_elast(pe,r) = 5;
$if "%prodspec%"=="resprod_06_gas15" supply_elast(pe,r) = 0.75 ; supply_elast("coa",r) = 0.8;  supply_elast("coa",r) = 1.5;

* formula from the appendix
nu(pe,r)$(((1 - beta_Res(pe,r)) and beta_all(pe,r)))= (supply_elast(pe,r) +1 - 1/beta_all(pe,r)) * (beta_all(pe,r) *  beta_Res(pe,r)    /  (1 - beta_Res(pe,r))) ;

res_sect(i)$(sum(r, beta_Res(i,r)))  = yes;

* for calculation of change in income

endow_sh_res(i,r)$res_sect(i) = vfm_adj("NatlRes",i,r) / sum((f.local,i.local),vfm_adj(f,i,r)) ;

CHECK_income(r) = (sum(f$f_noRes(f),endow_sh(f,r)) + sum(i$res_sect(i),  endow_sh_res(i,r) ));
* ok, sums to 1

*---------------------------------------------------------------------------------------------------

parameter bitc_shock;
* Trade cost shock just between different countries:
*bitc_shock(r,s) = tcost_shock;
bitc_shock(r,s) = 1;


positive variables
                chg_income(r) per capita income
                chg_lambda(r) lagrance multiplier
                chg_fd(i,r) final demand
                chg_absorption(i,r) absorption
                chg_production(i,r) production
                chg_Index(i,r) Price index P
                chg_supplier(i,r) Supplier effect S
                chg_fprice(f,r) factor price
*                chg_wage(f_Lab,r) wages

* only relevant for the sectors in res_sect
		chg_resprice(i,r)  price of sector specific resource
               
* intermediate nests to simplify notation (not necessary)
		res_nonres(i,r)   CES -- incl. income shock
		factorcost(i,r)	  C-D
		res_nonres_forfmc  sameas above but without income shock
 ;


equations eq_budget, eq_tradeprod, eq_index, eq_wage, eq_res_nonres(i,r), eq_factorcost , eq_FMC_nonres, eq_FMC_res, eq_res_nonres_forfmc ;

* Combining budget and demand:
$if "%spec%"=="NH" eq_budget(r)..       chg_income(r) =e= sum(i,dem_sh(i,r) * chg_lambda(r)**(-sigma(i)) * chg_index(i,r)**(1-sigma(i))  );

* if homothetic
$if "%spec%"=="H" eq_budget(r)..       chg_income(r) =e= sum(i,dem_sh(i,r) * chg_lambda(r)**(-1) * chg_index(i,r)**(1-sigma(i))  );


* Combining Trade, production and absorption:
eq_tradeprod(i,r)..  chg_production(i,r) =e=  sum(s,prod_sh_adj(i,r,s) * chg_supplier(i,r) * bitc_shock(r,s)**(-theta(i)) *
                         chg_index(i,s)**theta(i) * (gosh_f(i,s)*chg_fd(i,s) + sum(j$gosh_m(i,j,s),gosh_m(i,j,s)*chg_production(j,s))) );


* Combining index and supplier effect into a single equation:
eq_index(i,s)..      chg_index(i,s)**theta(i) * sum(r,trade_sh(i,r,s) * bitc_shock(r,s)**(-theta(i))  *
			prod(j,chg_index(j,r)**intersh(j,i,r))**(-theta(i)) *
			(res_nonres(i,r))**(-theta(i) * beta_all(i,r) ) 
			) =e= 1;

eq_res_nonres(i,r)..  res_nonres(i,r) =e=  (beta_Res(i,r) * ( (chg_resprice(i,r)* (1/income_shock))*income_shock_toresourcefactor
										+(chg_resprice(i,r))*(1- income_shock_toresourcefactor)
					 )**(1-nu(i,r)) + (1-beta_Res(i,r)) * (  factorcost(i,r) /income_shock   )**(1-nu(i,r)) )**(1/(1-nu(i,r)));


* income_shock is applied above
eq_factorcost(i,r) ..  factorcost(i,r) =e= prod(f$f_noRes(f),(chg_fprice(f,r))**beta_fnores(f,i,r));		
  

* cannot use the same nest for fmc's because we don't want the shock here..
* this is an intermdiate value to simplify notation
eq_res_nonres_forfmc(i,r)	 ..  res_nonres_forfmc(i,r) =e=  (beta_Res(i,r) * (chg_resprice(i,r) )**(1-nu(i,r)) + (1-beta_Res(i,r)) * (  factorcost(i,r) )**(1-nu(i,r)) )**(1/(1-nu(i,r)));



eq_FMC_nonres(f,r)$f_noRes(f) ..	chg_fprice(f,r) =e= sum(i$(not res_sect(i)),  fact_sh(f,i,r) * chg_production(i,r)) 
								+ sum(i$(res_sect(i)),  fact_sh(f,i,r) * chg_production(i,r)* 
								(factorcost(i,r)**(1-nu(i,r)))  * (res_nonres_forfmc(i,r)**(nu(i,r)-1))) ; 

eq_FMC_res(i,r)$res_sect(i) ..		chg_resprice(i,r)**nu(i,r) =e= chg_production(i,r) * (	res_nonres_forfmc(i,r)**(nu(i,r)-1));



* FOUR STEPS PER LOOP:
* 1) SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:
* 2) SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:
* 3) SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:
* 4) MANUALLY ADJUSTING FACTOR PRICES


* MODELS:

* 1) MODEL TO SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:
*model counterfactMCP_PE_justprices / eq_index.chg_index /;
model counterfactMCP_PE_justprices / eq_index.chg_index, eq_factorcost.factorcost, eq_res_nonres.res_nonres /;


counterfactMCP_PE_justprices.iterlim = 1000000;

* 2) MODEL SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:
model counterfactMCP_PE_demand / eq_budget.chg_lambda /;
counterfactMCP_PE_demand.iterlim = 1000000;

* 3) MODEL FOR SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:
model counterfactMCP_PE_fixedprice / eq_tradeprod.chg_production/;
counterfactMCP_PE_fixedprice.iterlim = 1000000;

* 4) MODEL FOR SOLVING FOR factor prices, TAKING PRODUCTION AS GIVEN:
model counterfactMCP_PE_wage / eq_FMC_nonres.chg_fprice, eq_FMC_res.chg_resprice,  eq_factorcost.factorcost, eq_res_nonres_forfmc.res_nonres_forfmc /;
counterfactMCP_PE_wage.iterlim = 1000000;


*------------

* APPROXIMATING STARTING VALUE FOR RESOURCE PRICE (speeds up convergence)



parameters  chg_fd_approx, chg_absorption_approx, chg_fprice_justfd, chg_fprice_approx, chg_skillprem_approx, chg_skillprem_justfd, log_Y_qty_approx;

chg_fd_approx(i,r) = 1 + (incelast("DIRECT",i,r) - 1) * (income_shock - 1);

display chg_fd_approx;


log_Y_qty_approx("approx",i,r ) =  log(income_shock) * (incelast("total",i,r))  ;
log_Y_qty_approx("approx",pe,r ) =  log(income_shock) * (incelast("total",pe,r)) * (1 - 1 /(1+ supply_elast(pe,r))) ;

display log_Y_qty_approx;


equation eq_FMC_res_approx;

eq_FMC_res_approx(i,r)$res_sect(i) ..		chg_resprice(i,r)**nu(i,r) =e= chg_production(i,r) * ( (beta_Res(i,r) * (
						chg_resprice(i,r)
						)**(1-nu(i,r)) + (1-beta_Res(i,r)) * (  factorcost(i,r)   )**(1-nu(i,r)) )**(-1));


model eq_FMC_res_presolve / eq_FMC_res_approx.chg_resprice /;

chg_production.FX(i,r) = ((income_shock-1) * supply_elast(i,r) + 1) /income_shock;



factorcost.FX(i,r) = 1;

res_nonres.L(i,r) =1;
chg_resprice.L(i,r) =1;


solve eq_FMC_res_presolve using MCP;

chg_production.UP(i,r) = 1000;
factorcost.UP(i,r) = 1000;

chg_production.LO(i,r) = 1/1000;
factorcost.LO(i,r) = 1/1000;

*------------------------


parameter count, Sqr_stat, Abs_stat;
count = 0;
Abs_stat = 100;
Sqr_stat = 100;

parameters      Sqr_stat_rep(T)
                Abs_stat_rep(T)
                chg_income_rep(T,r) per capita income
                chg_lambda_rep(T,r) lagrance multiplier
                chg_fd_rep(T,i,r) final demand
                chg_absorption_rep(T,i,r) absorption
                chg_production_rep(T,i,r) production
                chg_production_qty_rep(T,i,*) production
                chg_Index_rep(T,i,r) Price index P
                chg_supplier_rep(T,i,r) Supplier effect S
                chg_fprice_rep(T,f,r) factor price
                chg_resprice_rep(T,i,r) factor price
                chg_income_ini(r) per capita income
                chg_fprice_ini(f,r) factor price
                tradecosts_rep(r,s) trade costs
                chg_resprice_import(i,r) Storing resource factor price
                chg_fprice_import(f,r) Storing factor price;


tradecosts_rep(r,s) = 0;

chg_fprice_import(f,r) = 1;

chg_resprice_import(i,r) = chg_resprice.L(i,r);

bitc_shock(r,s) = tcost_shock;
bitc_shock(r,r) = 1;

tradecosts_rep(r,s)= bitc_shock(r,s);

display bitc_shock;

* Resetting count and convergence stats:
count = 0;
Abs_stat = 100;
Sqr_stat = 100;


*Bounds:
chg_income.LO(r) =  0.001;
chg_lambda.LO(r) =  0.001;
chg_fd.LO(i,r) =  0.001;
chg_absorption.LO(i,r) =  0.001;
chg_production.LO(i,r) =  0.001;
chg_Index.LO(i,r) =  0.001;
chg_supplier.LO(i,r) =  0.001;
chg_fprice.LO(f,r) =  0.001;

chg_resprice.LO(i,r) =  0.001;

res_nonres.LO(i,r) =  0.001;
res_nonres_forfmc.LO(i,r) =  0.001;


factorcost.LO(i,r) =  0.001;


chg_income.UP(r) = 1000;
chg_lambda.UP(r) = 1000;
chg_fd.UP(i,r) = 1000;
chg_absorption.UP(i,r) = 1000;
chg_production.UP(i,r) = 1000;
chg_Index.UP(i,r) = 1000;
chg_supplier.UP(i,r) = 1000;
chg_fprice.UP(f,r) = 1000;
chg_resprice.UP(i,r) = 1000;
res_nonres.UP(i,r) =  1000;
res_nonres_forfmc.UP(i,r) =  1000;
factorcost.UP(i,r) =  1000;

*Initial value:
chg_income.L(r) = 1;
chg_lambda.L(r) = 1;
chg_fd.L(i,r) = 1;
chg_absorption.L(i,r) = 1;
chg_production.L(i,r) = 1;
chg_Index.L(i,r) = 1;
chg_supplier.L(i,r) = 1;
chg_fprice.L(f,r) = 1;
chg_resprice.L(i,r) = 1;
res_nonres.L(i,r) =  1;
res_nonres_forfmc.L(i,r) =  1;

factorcost.L(i,r) =  1;


* INITIAL ADJUSTMENT FOR INCOME COUNTERFACTUAL:
* 
chg_Index.L(i,r) = 1/income_shock;
chg_fprice.FX(f,r) = 1;
chg_resprice.FX(i,r) = chg_resprice.L(i,r);
chg_Income.FX(r) = 1;

solve counterfactMCP_PE_justprices using MCP;

chg_index.FX(i,r) = chg_index.L(i,r);

chg_supplier.FX(i,r) =	 prod(j,chg_index.L(j,r)**intersh(j,i,r))**(-theta(i)) *
		(res_nonres.L(i,r))**(-theta(i) * beta_all(i,r) ) ;

display chg_index.L;

solve counterfactMCP_PE_demand using MCP ;

$if "%spec%"=="NH" chg_fd.FX(i,r) = chg_lambda.L(r)**(-sigma(i)) * chg_index.L(i,r)**(1-sigma(i));
$if "%spec%"=="H" chg_fd.FX(i,r) = chg_lambda.L(r)**(-1) * chg_index.L(i,r)**(1-sigma(i));

chg_lambda.FX(r) = chg_lambda.L(r);

display chg_lambda.L, chg_fd.L;

solve counterfactMCP_PE_fixedprice using MCP ;
display chg_production.L;

chg_absorption.L(i,s) = gosh_f(i,s) *chg_fd.L(i,s) + sum(j$gosh_m(i,j,s),gosh_m(i,j,s) *chg_production.L(j,s));

* T LOOPS:
loop(T,

*for counterfactual in tau:
If ((Abs_stat > 0.0000001),

* 1) SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:

chg_supplier.LO(i,r) = 0.001;
chg_supplier.UP(i,r) = 1000;
chg_index.LO(i,r) = 0.001;
chg_index.UP(i,r) = 1000;

solve counterfactMCP_PE_justprices using MCP ;


chg_index.FX(i,r) = chg_index.L(i,r);

*chg_supplier.FX(i,r) = prod(f$(f_noLab(f) and beta(f,i,r)),(chg_fprice.L(f,r)/income_shock)**beta(f,i,r))**(-theta(i)) *
*                       sum(f$(f_Lab(f) and beta(f,i,r)),(beta(f,i,r)/betalab(i,r))*(chg_fprice.L(f,r)/income_shock)**(1-rho))**(-theta(i)*betalab(i,r)/(1-rho)) *
*                       prod(j,chg_index.L(j,r)**intersh(j,i,r))**(-theta(i));

chg_supplier.FX(i,r) =	 prod(j,chg_index.L(j,r)**intersh(j,i,r))**(-theta(i)) *
			(res_nonres.L(i,r))**(-theta(i) * beta_all(i,r) ) ;


* 2) SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:

chg_fd.LO(i,r) = 0.001;
chg_fd.UP(i,r) = 1000;
chg_lambda.LO(r) = 0.001;
chg_lambda.UP(r) = 1000;

solve counterfactMCP_PE_demand using MCP ;

$if "%spec%"=="NH" chg_fd.FX(i,r) = chg_lambda.L(r)**(-sigma(i)) * chg_index.L(i,r)**(1-sigma(i));
$if "%spec%"=="H" chg_fd.FX(i,r) = chg_lambda.L(r)**(-1) * chg_index.L(i,r)**(1-sigma(i));

chg_lambda.FX(r) = chg_lambda.L(r);

* 3) SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:

solve counterfactMCP_PE_fixedprice using MCP ;

chg_absorption.L(i,s) = gosh_f(i,s) *chg_fd.L(i,s) + sum(j$gosh_m(i,j,s),gosh_m(i,j,s) *chg_production.L(j,s));


* 4) SOLVING FOR WAGES, TAKING PRODUCTION AS GIVEN:

chg_production.FX(i,r) = chg_production.L(i,r);

chg_fprice.LO(f,r)$f_noRes(f) = 0.001;
chg_resprice.LO(i,r)$res_sect(i) = 0.001;

chg_fprice.UP(f,r)$f_noRes(f) = 1000;
chg_resprice.UP(i,r)$res_sect(i) = 1000;


solve counterfactMCP_PE_wage using MCP ;

chg_production.LO(i,r) = 0.001;
chg_production.UP(i,r) = 1000;

* 4) CONVERGENCE STATS + MANUALLY ADJUSTING FACTOR PRICES:

Sqr_stat = ( sum((r,f)$f_noRes(f), sqr(1-  chg_fprice.L(f,r) /  chg_fprice_import(f,r)   )) /  sum((r,f)$f_noRes(f),1)
		+ sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)), sqr(1-  chg_resprice.L(i,r) /  chg_resprice_import(i,r) ) ) /  sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)),1) )**(1/2) ;


Abs_stat = sum((r,f)$f_noRes(f), abs(1-  chg_fprice.L(f,r) /  chg_fprice_import(f,r) ) ) /  sum((r,f)$f_noRes(f),1) 
		+ sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)), abs(1-  chg_resprice.L(i,r) /  chg_resprice_import(i,r) ) ) /  sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)),1)  ;


count = count + 1;

display Sqr_stat, Abs_stat, count, chg_fprice.L;


* weights can be 1/3 and 2/3rd for income counterfactual.
*chg_fprice.FX(f,r) = (1/3) * chg_fprice.L(f,r) + (2/3) * sum(i,fact_sh(f,i,r) * chg_production.L(i,r));


chg_fprice.FX(f,r) = (1/5) * chg_fprice.L(f,r) + (4/5) * chg_fprice_import(f,r);


* 1/10 works well  (2/10 doesnt work)
chg_resprice.FX(i,r) = (1.5/10) * chg_resprice.L(i,r) + (8.5/10) * chg_resprice_import(i,r);

* compute change in income here

chg_income.FX(r) = sum(f$f_noRes(f),endow_sh(f,r)*chg_fprice.L(f,r))
	 + sum(i$res_sect(i),  endow_sh_res(i,r) * chg_resprice.L(i,r));


chg_fprice_import(f,r) = chg_fprice.l(f,r);
chg_resprice_import(i,r) = chg_resprice.l(i,r);

	

* Reporting:
Sqr_stat_rep(T) = Sqr_stat;
Abs_stat_rep(T) =  Abs_stat;
chg_lambda_rep(T,r) = chg_lambda.L(r);
chg_fd_rep(T,i,r) = chg_fd.L(i,r);
chg_absorption_rep(T,i,r) = chg_absorption.L(i,r);
chg_production_rep(T,i,r) = chg_production.L(i,r);
chg_Index_rep(T,i,r) = chg_Index.L(i,r);
chg_supplier_rep(T,i,r) = chg_supplier.L(i,r);
chg_fprice_rep(T,f,r)$f_noRes(f) = chg_fprice.L(f,r);
chg_resprice_rep(T,i,r)$res_sect(i) = chg_resprice.L(i,r);
chg_income_rep(T,r) = chg_income.L(r) ;
chg_production_qty_rep(T,i,r) = chg_production_rep(T,i,r) / chg_Index_rep(T,i,r);
chg_production_qty_rep(T,i,"world") = sum(r, production(i,r) * chg_production_rep(T,i,r) / chg_Index_rep(T,i,r)) / sum(r, production(i,r));

);
* end of "if" statement

display count;

);
* end of loop across T's


* dump all parameters 
execute_unload 'results/CFresults_co2_%spec%_%demandest%_%prodspec%_fitteddem.gdx'
*execute_unload 'results/CFresults_co2_%spec%_%demandest%_%prodspec%_fitteddem_10pctshock.gdx'


